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Understanding the non-equilibrium dynamics of extended quantum systems after the trigger of 
a sudden, global perturbation (quench) represents a daunting challenge, especially in the presence 
of interactions. The main difficulties stem from both the vanishing time scale of the quench event, 
which can thus create arbitrarily high energy modes, and its non-local nature, which curtails the 
utility of local excitation bases. We here show that nonperturbative methods based on integrability 
can prove sufficiently powerful to completely characterize quantum quenches: we illustrate this 
using a model of fermions with pairing interactions (Richardson's model). The effects of simple 
(and multiple) quenches on the dynamics of various important observables are discussed. Many of 
the features we find are expected to be universal to all kinds of quench situations in atomic physics 
and condensed matter. 



The experimental realization [lj of cold atomic systems 
with a high degree of tunability of Hamiltonian param- 
eters, and the ability to evolve in time with negligible 
dissipation, has reignited the study of many-body quan- 
tum systems away from equilibrium. How Gibbs or any 
other relaxed states can ultimately result from unitary 
dynamics is a question that has received a lot of atten- 
tion recently [2J |3J SJ EJ E] , but which still lacks a general 
understanding. 

Suppose an extended quantum system is prepared in 
one eigenstate \i/jg ) of some Hamiltonian H go , where go 
is a tunable, global parameter (interaction strength, ex- 
ternal field, ...). At a given time, say t = 0, this pa- 
rameter is suddenly changed to a different value g, and 
the system thus starts evolving unitarily according to the 
dynamics governed by a different Hamiltonian H g . This 
is what is referred to as a quantum quench. The result- 
ing time evolution is simply given by the solution of the 
Schrodinger equation \ip(t)) — e~ iHat \ipg )- Since li^gg) 
is not an eigenstate of H g , this can be extremely diffi- 
cult to quantify. The most straightforward way to tackle 
the problem is therefore to write the initial state \ipg ) 
as a sum over the complete set of eigenstates \i/jg) (hav- 
ing energy ujg) of Hamiltonian H g , leading to the time- 
dependent post-quench state 

V 

The complexity of the problem is encoded first in the 
distribution of energies uj v g , but most importantly in the 
matrix of overlaps of eigenstates pertaining to the two 
different Hamiltonians, 

0%, = > ( 2 ) 

which we call the quench matrix. This matrix is of di- 
mension equal to that of the Hilbert space [51] . but in 



practice we mainly need a single column expressing the 
initial eigenstate of H go in terms of the eigenstates of H g . 
However, even in the very few cases when all eigenstates 
of a many-body Hamiltonian can be classified and written 
down, the calculation of the quench matrix coefficients is 
a severe challenge, whose computational complexity gen- 
erally grows factorially with system size. Shortcuts can 
be found for systems having a representation in terms 
of free particles (like the ID Ising chains [5j [6]) where 
Wick's theorem suffices to calculate all the overlaps, but 
for truly interacting systems this remains a very ambi- 
tious programme. Most of the theoretical work on quan- 
tum quenches has up to now concentrated on the calcu- 
lation of correlation functions in specific regimes [2 [3] , 
with little reference to the post-quench state of the sys- 
tem. 

Besides describing the state resulting from a quench of 
state /i, it is also important to be able to characterize 
the time dependence of physical observables O after the 
quench. Formally, we can write 

(o(t)) = mm*®) 

= E^-^W^o^lOl^), (3) 

where calculating the matrix elements (ipg] O \ipg) rep- 
resents an additional hurdle for interesting observables 
in nontrivially interacting models. Even if we are able 
to obtain these matrix elements, the leftover double sum 
over the full Hilbert space is enormous, and one can won- 
der whether this way of proceeding is of any practical use. 
New nonperturbative methods are clearly needed to ob- 
tain a proper description of the physics involved. 

The purpose of the present paper is to introduce a new 
line of attack on quantum quench problems, sufficiently 
powerful to yield not only the quench matrix of specific 
interacting problems (and thus the ensuing nonequilib- 
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rium state), but also able to provide matrix elements of 
physical observables, and thus their time dependence af- 
ter the quench. This approach is based on the exact solv- 
ability of certain many-body quantum problems known 
as integrable or Bethe Ansatz [7J solvable theories. Inte- 
grability came into prominence as a means of obtaining 
exact results for the equilibrium thermodynamics of one- 
dimensional systems (see [8j [9J and references therein). 
More recently, a description of correlation functions at 
equilibrium has been achieved by exploiting results from 
the Algebraic Bethe Ansatz (ABA), which provides eco- 
nomical expressions for matrix elements of local opera- 
tors in the basis of exact Bethe eigenstates. The existence 
of these expressions stems from two results: Slavnov's 
formula [TU] for the overlap of a Bethe state with a generic 
state, and the solution of the so-called quantum inverse 
problem i.e. the mapping of physical operators to 
ABA operators. These matrix elements are of great util- 
ity in the computation of equilibrium correlation func- 
tions. One very important feature is that Bethe states 
typically offer a very optimized basis in which only a very 
small minority of eigenstates carry substantial correla- 
tion weight, allowing the summation over intermediate 
states to be drastically truncated without significantly 
affecting the results. This novel approach has been suc- 
cessfully applied to equilibrium correlations of quantum 
spin chains [T3] and atomic Bose gases [T2]. We here fur- 
ther extend the reach of integrability into the domain of 
non-equilibrium quench dynamics. 



THE MODEL AND ITS SOLUTION 

We consider a model of spin 1 /2 fermions in a shell of 
energy levels e a with a Cooper pairing-like interaction 

ff = E I! y c «' C[w ■ j I! c i+ c i- c f- c w ■ ( 4 ) 

a a a, (3 

which was introduced by Richardson |14j in the context of 
nuclear physics, and has found applications in the physics 
of ultrasmall metallic grains 15j. It reduces to conven- 
tional Bardeen-Cooper-Schrieffer (BCS) theory [TO] in 
the thermodynamic (TD) limit. The model has a pseudo- 
spin representation with S~ — c Q |C Q | (see Methods) with 
N spins. Central to our approach is the fact that this 
Hamiltonian can be diagonalized using the Bethe Ansatz 
[TU [17] . The Hilbert space separates into sectors of fixed 
number of down spins N r . The eigenstates of the model 
are given by Bethe wavefunctions, each individually char- 
acterized by a set of N r rapidities {wj} obeying a set of 
algebraic equations known as the Richardson equations 

N N r 

l = J2^— -E^— i = l>-.*r, (5) 
q i — ' Wj — e„ ' — ' Wj — Wi? 

y Q= l J k^j 1 K 



and are obtained by repeated action of an operator B(wj) 
on the fully polarized reference state |0): 

2V r N 

ik}) = n*K-)io> = n e — m°>- ^ 

J k — 1 a—1 

The ) different solutions to (p| then allow to construct 
a full set of orthogonal eigenstates, providing us with a 
proper basis of the Hilbert space. 

In Bethe Ansatz solvable models, Slavnov's formula 
[TO], gives the overlap of an eigenstate \{w}) with a gen- 
eral Bethe state \{v}) (see Methods). The main difference 
with other models solvable by Algebraic Bethe Ansatz 
(like the one-dimensional Bose gas or the XXZ chain) 
is that in the definition of the eigenstates pj, the cou- 
pling constant g enters only implicitly through the solu- 
tions of the Richardson equation for Wj . Consequently, 
for this model, Slavnov's formula is enough to calculate 
the overlaps between two generic states at any coupling. 
For other models, where B{wj) depends explicitely on g, 
it can only be used between states defined by the same 
operators B(wj) and therefore the same g, and a more 
general expression for the overlaps remains to be found. 

We can then exploit the accessibility to the quench ma- 
trix for the Richardson model to show how useful integra- 
bility can be when studying quenches. However, before 
entering in the details of the quantum dynamics, it is im- 
portant to remember that in the TD limit N r ,N — ► oo 
at fixed filling, the dynamics becomes classical [HI [20] 
because of suppression of quantum fluctuations. In this 
limit, the dynamics of the canonical order parameter can 
be obtained analytically by exploiting classical integra- 
bility [TO] [50]. The framework we propose in this let- 
ter works in the mesoscopic regime (finite N), allowing 
to study the effects of quantum fluctuations. With this 
tool at hand, we can characterize in an exact manner the 
crossover taking place between microscopic and macro- 
scopic physics, a task impossible to achieve with thermo- 
dynamics! approaches. 

We will argue in the following that, similarly to what 
is observed for equilibrium correlation functions, only a 
relatively small set of states contributes significantly to 
the decomposition of the initial state in the new eigen- 
basis. The natural approach is then to truncate in an 
optimal way the Hilbert space, so that a faithful repre- 
sentation of the initial state is obtained. The induced 
truncation error is easily evaluated looking at how close 
KV'gJV'g)! is to the desired value of 1. Using the 
truncated Hilbert space we can then calculate any ob- 
servable or correlation function by brute force summing 
the relevant contributions. 

Compared to other numerical truncation methods, this 
approach has the great advantage that time enters only 
as a parameter. The explicit expression (JT|) for the wave- 
function means that at any time, expectation values can 
be computed without knowing the previous history of the 
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FIG. 1: First column of the quench matrix (ground-state over- 
laps) for several quenches. In all plots N = 2N r — 16 and the 
ground state energies (represented by vertical lines) have been 
shifted for clarity. Top: Decomposition of the g — ground- 
state with states at g — 0.05,0.5,0.95. Center: Decompo- 
sition of several initial ground-state go = 0,0.15,0.3,0.5 in 
terms of the states at g = 1. Bottom: Decomposition of the 
go = 1 ground-state in terms of g = 0.95, 0.55, 0.15, states. 



system (apart from the initial state) and there is therefore 
no accumulation of errors as time passes. On the other 
hand, compared to numerical exact diagonalization, ma- 
trix elements can be expressed using Slavnov's formula 
as matrix determinants whose computational complex- 
ity is algebraic and not exponential in the system size 
and/or number of excitations, thereby allowing to reach 
large system sizes. 



NUMERICAL RESULTS 

Let us now report explicit results starting with the 
quench matrix itself. We concentrate on the case of 
equally spaced levels e a at half-filling, but the method 
is clearly not limited to this case. Let us start from 
TV = 2N r = 16, when the Hilbert space has a dimen- 
sion of "only" 12870. We can then numerically compute 
the rapidities for every single state. We report in Fig. 
[I] the square of the overlaps for several quenches as ob- 
tained from Slavnov's formula (see Methods). Starting 
from the non-interacting (g = 0) ground state, the top 
inset shows the overlaps with all the states at three dif- 
ferent finite couplings. This allows to understand some 
general features: having access to the complete quench 
matrix, one realizes that only a few of the eigenstates 
at coupling g have a large contribution to \ipo). There- 
fore, getting a nearly exact description of the dynamics 
requires only a small subset of the states. The final state 
having the largest overlap with the initial ground state 
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FIG. 2: Left: Pictorial representation of single block states 
obtained by promoting contiguous blocks of N p < N r rapidi- 
ties from right below to right above the Fermi level. Right: 
Total contribution of these states to the amplitude of the ini- 
tial state at on = 0, as a function of interaction. 



is always one of the states that at g = is built by flip- 
ping from up (down) to down (up) the N p spins right 
below (above) the Fermi level. We refer to these N r + 1 
states as 'single block states'. For the cases we stud- 
ied here, quenching from a non-interacting initial state, 
these states always contribute more than 60% of the total 
amplitude (see Fig. [2). The total contribution of single 
block states is non-monotonic in g showing (after decline 
at small interaction) a rise as g gets sufficiently large. 

The 'band-like' structure of the overlaps in Fig. [T] 
makes it reasonable to assume that additional large con- 
tributions can be found for states built by slightly de- 
forming the single block ones, e.g. by adding a single 
particle-hole excitation cither above or below the Fermi 
level. The remaining most relevant states will then be 
those with two additional excitations. Following this as- 
sumption inductively, we add to the truncated Hilbert 
space multiple block states obtained by slightly deform- 
ing single block ones. This allows us to describe larger 
systems while retaining a tractable number of states. 
This procedure works extremely well, e.g. at N = 32, 
for all the quenches from go = to g € (0, 1], we were 
always able to find at minimum 97% of the weight of the 
initial state by using only 7000 states, i.e. only 1/10 5 of 
the full Hilbert space. Moreover, for a given final value 
of g, less than a 1000 of these states gives an actual im- 
portant contribution. 

In the center of Fig. [T] we report the overlaps obtained 
by quenching from different initial values of go to the 
same final g = 1. We see that the same band- like struc- 
ture is present as when starting from a non-interacting 
state, leading to no qualitative change of the dynamics. 
Vice-versa, the structure of the quench matrix for a re- 
versed quench, i.e. from large to small <?, reported in the 
bottom of Fig. [T]is different: the weight of the state goes 
down exponentially with the energy of the states (almost 
straight lines in the figure), and for large quenches the 
decay rate in energy is slow. An adequate representation 
of the initial state therefore requires that lots of states 
be taken into account. In this case, the optimal trunca- 
tion procedure is still easily defined by simply keeping a 
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of the state with the right quantum numbers at the given 
energy. In the top panel we report several quenches at 
N = 32, where the formation of subdominant peaks is 
explictly shown. In the bottom panel we show P(W) 
at different N keeping gN fixed. It is evident that de- 
spite the fact that the structure changes drastically with 
N, the width of the distribution is constant, indicating 
that, when written in terms of W/N, it becomes a delta 
function. 
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FIG. 3: Distribution of work P(W). Top: Different quenches 
at fixed N = 32, showing the formation of multiple peaks for 
large quenches. Bottom: At different N keeping fixed gN. 



sufficient number of low energy states. 

The first measurable quantity easy to derive from the 
knowledge of the quench matrix is the probability distri- 
bution of the 'work' [B] 



(7) 



reported for some quenches from g = in Fig. [3] These 
figures have been obtained by smoothing the 5 function 
with a Gaussian of width proportional to the inter-level 
spacing. It is straightforward to derive the average and 
the width of this distribution 
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diagonal order-parameter in the initial state, and ^ 2 
(V'soKEa./^i S aSp/N r ) 2 \iP° go ) is a four point correlator. 
Higher cumulants are similarly obtained and only de- 
pend on the initial state [25]. From TD relations [6], it 
is generically expected that the probability of work per 
spin w = W/N becomes a delta function. For finite N, 
P(W) is non-trivial: it shows a dominant peak close to 
W = (W), but with a structure dictated by the presence 



We now present results for observables, concentrating 
on the off diagonal order parameter defined as 



*oD(t) = m)\ 



N, 



N 

E 

a,/3=i 



(9) 



In the equilibrium canonical ensemble ton for N — > oo 
is related to the BCS gap and so it is a natural quan- 
tity to understand the superconducting tendency even 
out of equilibrium (on the same footing as the canoni- 
cal order parameter used in [H]). According to Eq. ^ 
we can write the time evolution once the form factors 
for J^a /9=i &a are known. They have a representa- 
tion in terms of a sum of N r determinants of N r x N r 
matrices depending on the rapidities (see Methods). In 
the bottom left part of Fig. [4] we present the result- 
ing real-time evolution of ^oo(t) starting from g = 
and evolving with several different g for TV = 32. The 
information contained here is better extracted from the 
Fourier transforms reported on the right of Fig. |4j 

For small value of g, the various frequencies entering 
are very close to integers, as a result of the almost per- 
fect equispacing of the levels. This regime could simply 
be described by perturbation theory and doesn't show 
any striking features different from free fermions. With 
increasing g, the spectrum becomes very complicated 
since a large number of incommensurate frequencies con- 
tribute to the order parameter evolution. This is the 
realm of quantum fluctuations which makes the evolu- 
tion highly irregular. Still increasing g, some regularity 
appears again. This can be understood in terms of re- 
sults in the TD limit [T^]. In fact, for N — > oo, quench- 
ing from weak interaction to a much larger one, leads to 
an order parameter which shows persistent non-harmonic 
periodic evolution, i.e. a Fourier transform with equis- 
paced peaks. Within the canonical description presented 
here, this feature will be reproduced when quenching to a 
large final value of g since, as was shown in [5T], excited 
states in this regime form equispaced bands at energy 
E ps AN G + O(N ) (N G being the number of Gaudinos, 
i.e. the relevant excitations in this regime [H]). As can 
be readily seen looking at the energy distribution of the 
points in Fig. [T] for finite large couplings, the low energy 
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FIG. 4: Bottom left: Off-diagonal order parameter evolu- 
tion for N = 32. Right: Fourier transform, the various plots 
are shifted on the vertical axis for clarity. Top Left: Non- 
equilibrium finite-size "phase diagram" resulting from the 
time-averaged canonical gap obtained from the off-diagonal 
order parameter, as explained in the text. 



bands are already clearly formed, progressively collaps- 
ing into a single energy E « ANq. The slight remaining 
width of these bands would only result in additional low 
frequency corrections to the mean- field BCS result. 

The static correlation functions studied in [18 , 22] de- 
pend mainly on low energy properties, and BCS- like be- 
havior was always found when g > g* = (21niV) _1 , i.e. 
the criterion for the presence of superconductivity. In the 
problem at hand though, the quench matrix clearly shows 
that quantum quenches lead to an important occupation 
of the higher energy bands. These clearly differ from the 
BCS spectrum even for values of g much larger than g* 
(that for N — 32 is only 0.144 . . . ). Quenches, since they 
probe high energy properties of the system, open up the 
possibility of probing interaction effects not captured by 
the mean- field treatment. Non mean-field features arc 
manifest in non-equispaced dominant peaks in the fre- 
quency dependence of the order parameter. These ac- 
cessible experimental quantities could thus be used as a 
spectroscopic tool (as proposed for other models in Ref. 
[4J to study quantum fluctuations. 

We also considered the evolution of the canonica l order 
parameter defined as ^(t) = J2a=i \J\~ i^oSf)) 2 ^ us ing 
the knowledge of the form factors for S* (see Methods) . 
It displays the same qualitative features as *&oD(t) and 
consequently will not be discussed here. 

Let us conclude this section with a discussion of the 
long time asymptotic. It is difficult to extract any infor- 
mation about it from the highly irregular and oscillatory 
behavior reported in Fig. [4] Furthermore in finite sys- 
tem, (approximate) quantum recurrence will always spoil 
any signature of an eventual asymptotic state. However, 
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FIG. 5: Left: A typical quench go — > g populates all states 
according to the quench matrix. A general dequench g — > go 
redistributes the weight among all states. Right: A targeted 
dequench after a chosen time can populate a targeted state. 



if in the TD limit the asymptotic value of an observ- 
able exists, it must be equal to its time average, that 
is straightforwardly obtained with the tools at hand in 
finite systems. In fact, in Eq. ^ all the terms with 
v ± i average to zero and so O = £„ \Q g » g \ 2 ^ g \O\r g ), 
(the over line stands for the time average). $od is ob- 
tained with little effort using the Hellmann-Feynman the- 

orem (V g \ E^=i S+Sj \%) = -^S thus without in- 
volving determinants. The non-equilibrium 'phase dia- 
gram' for N — » oo [19 shows that the final asymptotic 
value of the canonical gap Aqo defines a universal curve 



where A g is the 



when expressing A^/Ag vs A go /A g 
equilibrium value at coupling g. In our normalization, 
A g = g^/Nr = g^/^oD/Nr - l/N r [18 (the -1 cancels 
the first correction for large N r ). We take the last equa- 
tion also off-equilibrium for the definition of the time av- 
eraged canonical gap in finite size. The resulting 'finite- 
size phase diagram' is reported in Fig. [4] (top left), where 
most of quenches from g to g both in [0, 1] are shown, 
for N — 8, 16, 32 at half filling (we excluded the points 
with an equilibrium ^od much different from BCS pre- 
diction, that are not expected to approach the asymptotic 
result). It is evident how increasing N the curves tend 
to the Barankov-Levitov result shown as a full line 261. 



THE DOUBLE QUENCH 

We move now to address the very interesting dynam- 
ics that appear when we consider sequences of multiple 
quenches, gi ~ > gi+\ at times ti. For brevity, we con- 
centrate here on the problem of the double quench, or 
quench-dequench sequence, defined by g — g for t < 0, 
g = gi for < t < t q and g — go for t > t q . Starting 
from a specific eigenstate of H go , the quench at t = Q 
populates excited states of H gi according to the quench 
matrix ([2]). Letting the system evolve up to time t q and 
then 'dequenching' back to go results in a nontrivial am- 
plitude of occupation for eigenstates of H go , given by the 
quench propagator 



P/3a (tq 



E 



Qgogi ( 



■gigo 



(10) 



6 



where a, ft S H go are respectively the labels for the pre- 
and post-quench states. For t q = 0, this propagator falls 
back onto the identity matrix. For finite duration t q > 0, 
interference effects lead to nontrivial states (see left panel 
of Fig. [5]). For a specific initial state a and final state /3, 
the quench propagator can be visualized as the sum of 
arrows of length | Q@ ^ Q2™ go | rotating as a function of t q 
at frequency ui^ from an initial phase arg(Q^ 9i Q^ i Q ffo ). 
When arrows of non-negligible length align, constructive 
interference occurs, favoring the weight of the final state 
(3 (see right panel of Fig. [5]). Since all arrows rotate at 
different frequency, the occupation probability of state (3 
is a highly nontrivial function of the quench time, which 
is however completely characterized from the information 
we now have at hand. 

We consider for dcfinitcncss a double quench start- 
ing from the ground state of Hamiltonian H go . As a 
function of the quench duration t q , the amplitudes of 
eigenstates a after the dequench will thus be given by 
A a {t q ) — P a o(t q ). We present in Fig. [6] the results of 
such double quench calculations. We specifically use a 
system of 16 spins, and trace over all intermediate states, 
allowing us to verify that the sum of square amplitudes 
remains equal to one (up to numerical accuracy of or- 
der 10 -7 ) at all quench times. The top panel shows the 
ground state occupation, which is inevitably the domi- 
nant state for small t q . However, we surprisingly find 
that its weight essentially vanishes (square amplitude be- 
low 0.005), first around quench time t q — 0.56, and also 
repeatedly afterwards. The ground state is also periodi- 
cally reconstructed to a large degree, showing that sub- 
stantial sloshing of the occupation weight in the Hilbert 
space occurs as a function of the quench duration. 

The occupation of individual excited states after the 
double quench also displays prominent time-dependent 
interference effects. Their amplitudes all begin at zero 
for t q = 0, but individual states can attain non-negligible 
amplitudes when the 'arrows' in their quench propagator 
add up constructively for particular quench durations. 
This is shown in the middle panel of Fig. |6j where we 
plot the occupation probability of three example states 
among the single block states. The times at which such 
alignments take place can be predicted using a simple 
algorithm based on what could be called a continuous 
sieve of Eratosthenes. Namely, for a given final state /3, 
the double quench weights \Q g g 9l QJ" go \ f° r au 7 are nrs t 
ordered in decreasing value. The dominant mode (rela- 
beled 0) has a time-dependent phase 4>°{t q ) = u gi t q — </>° 

with (f>° = aig(Q^g i Q g \g o ) , with similar defined phases 
for the subdominant modes i > 0. Choosing an arbi- 
trary phase alignment tolerance 89, the requirement that 
\cj) l (t q ) — (j) (t q )\ < 89 for agiven 'arrow' i defines excluded 
time intervals on the quench timeline t q € [0, 00 [. Erasing 
all such intervals for all states up to a level n leaves only 
the times at which all phases (j)°(t q ), 4> n (t q ) are aligned 
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FIG. 6: Occupation probabilities and moments after a double 
quench between g = and 0.5, as a function of the quench 
duration t q , for a system of 16 spins. The top plot gives the 
ground state occupation probability, clearly displaying van- 
ishing and reconstruction effects. The middle plot gives the 
occupation probabilities for the three more relevant states. 
The lower graph shows the moments I q , quantifying the de- 
gree of localization in Hilbert space after the double quench. 

to the chosen tolerance, and for which a certain amount 
of constructive interference occurs. For example, in the 
middle panel of Fig. |6j the b state peaks around t q ~ 4.8; 
it can be checked that this is a level 8 alignment (with 
tolerance chosen as 89 = 7r/8). Alignments of a given 
order n and tolerance 89 occur more or less periodically. 
Increasing the order or reducing the tolerance 89 makes 
alignments of higher quality but quickly increasing rarity. 
In view of this sieve of Eratosthenes logic, an interesting 
question is whether the distribution of quench alignment 
times can be linked to that of e.g. prime numbers. 

A study of the amplitudes A a after a double quench for 
each individual final state is clearly prohibitive. Charac- 
terizing the distribution of amplitudes is more enlight- 
ening, and is best performed by exploiting tools com- 
mon in the theory of localization in disordered systems, 
i.e. by considering the inverse participation ratios (IPRs) 
I q = J2 a \ A a \ 2q , with 1\ = 1. In the bottom panel of Fig. 
[6] we plot the second and third IPRs for excited states 
(defined as V = £ Q>0 K| 2 7(E Q >o \A a \ 2 )^ i.e. sum- 
ming over excited states only), which display the localiza- 
tion tendencies of the excited states' amplitude weight in 
the Hilbert space after the double quench. Curves of l2, r 
and Iz t r (always < l2,r) approach one another when one 
excited state becomes dominant, and indicate smoother 
weight distribution otherwise. 

Another interesting quantity to look at, which also has 
the advantage of being more directly accessible in exper- 
iments, is the work 

w{t q ) = Y,«-0\Mt q )\ 2 (ii) 
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FIG. 7: Total energy (work) pumped into the system by the 
quench-dequench sequence with go — and gi = 0.5 and 16 
spins, as a function of the quench duration t q . Inset: work as 
a function of t q for the first oscillations of the envelope. These 
persist for much longer times, which are not plotted for clarity. 
Main plot: Fourier transform of the work, showing the main 
peak associated to the energy difference of the two dominant 
intermediate states. 

or in other words the energy which is pumped into the 
system by a quench-dequench sequence of duration t q . 
Starting from the ground state means that W(t q ) is 
strictly positive. Since the quench-dequench sequence 
populates excited states in a highly t q dependent way, 
this quantity will also display a rich frequency profile. 
In Fig. [7] we plot (inset) the work as a function of t q , 
which displays nontrivial oscillatory behavior dominated 
by a frequency oj ~ 4.62 corresponding to the energy 
difference between the two dominant intermediate states 
during the quench. The Fourier transform W(u q ) is plot- 
ted in the main part of the figure, clearly displaying the 
above-mentioned peak but also the non-negligible contri- 
butions from a broad range of frequencies. The position 
of the peaks corresponds to excited energy level differ- 
ences of the Hamiltonian during the quench, their height 
giving information on the size of the relevant quench ma- 
trix elements. The work can thus be used not only as a 
spectroscopic tool, but as a way to quantify eigenstate 
overlaps. 



DISCUSSIONS 

In this work, we have proposed a novel method to 
tackle quantum quenches, based on integrability. We ap- 
plied the method to the fermionic pairing model, showing 
the large amount of information that can be obtained on 
e.g. the work probability density, physical observables 
(and correlation functions) and their time evolution, and 
multiple quenches settings. Everything has been derived 
in an exact (or numerically exact) manner in the meso- 
scopic regime, where quantum fluctuations govern the 
dynamics. We gave evidence of how single and double 
quench dynamics can be effectively used to extract spec- 
troscopic data from simple measurable quantities like the 



work done on the system. 

To obtain these results we explored the peculiar prop- 
erty that the quench matrix of the pairing model can be 
obtained using Slavnov's formula. This is not true in a 
general integrable model, and the generalization of the 
quench matrix representation is an open problem in the 
theory of integrable systems. When this representation 
will be available, the methods we propose here will al- 
low exact calculations for a large variety of experimental 
relevant models, most importantly the one-dimensional 
Bose gas and Heisenberg spin chains. 
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METHODS 
Solving the Richardson equations 

Because of the blocking effect excluding singly occu- 
pied levels from the dynamics [14], Richardson's model 
also has a pseudospin representation S~ — c a ic a ^, 5+ = 

c aT c ai' $a = C L| C al Ca i Cq T ~ V 2 - Tne Hamiltonian be- 
comes 

N N 
a=l a, 0=1 

where N is the number of unblocked levels. 

At g = the (^) solutions to the Richardson equa- 
tions are trivial. They are given by Eq. ([6| with the 
N r rapidities set to be strictly equal to one of the ener- 
gies e a . Apart from a few particular cases with a small 
number of particles, the Richardson equations are not 
solvable analytically when g ^ 0, and one should solve 
them numerically. The solutions are such that every Wj 
is either a real quantity or forms, with another parameter 
Wj>, a complex conjugate pair (CCP), i.e. w*, = Wj. The 
mechanism for the CCPs formation is very easy: as inter- 
actions are turned on, all Wj are real quantities for small 
enough g, but at a certain critical value of the coupling g* 
two rapidities will be exactly equal to one of the energy 
levels (wj = Wj/ = e 7 (j)) and for g > g*, the two param- 
eters that collapsed will form a CCP at least for a finite 
interval in g. The situation is in fact rather intricate: 
the values g* are implicit functions of all other rapidities, 
and can only be read off a full solution of the Richard- 
son equations for a specific choice of state. Moreover, 
CCPs can split back into real pairs, whose components 
can then re-pair with neighboring rapidities. Different 
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choices of the parameters e Q and of their eventual degen- 
erations specify different models. We specialize to the 
case of equally spaced levels. We make the choice to use 
e a = a which sets the zero of energy and implies that ev- 
ery energy will be given in units of the (pair) inter-level 
spacing. Furthermore we consider only half-filling of the 
energy levels (N = 2N r ), when the number of rapidities 
N r equals the number of particles N p , while in general 
in our notations N r + N p — N. At the precise value of 
g at which a pair of rapidities (wj,Wji) collapse into a 



CCP { Wj 



w 



e 7 (j')), the Richardson equations (|5|) 



labelled j and f will include two diverging terms whose 
sum remains finite. In order to be able to treat these 
points numerically, one can define the real variables, 
Wi j = Wj+Wji and W2,j = (2e 1 (j)—Wj—Wj>)/(Wj—Wj>). 
As first discussed in [53], we need to know beforehand 
which rapidities will form a CCP and at which e 7 (j) it 
will happen in order to use this type of change of vari- 
ables. Here we find the various solutions to the Richard- 
son equations numerically by increasing g by small steps 
starting from the solution at g = and can therefore pre- 
dict at every step, the upcoming formation of CCPs. As 
a consequence of this procedure, any given state at finite 
coupling can then be defined uniquely by the g — state 
from which it emerges (and the actual value of g). 



Scalar products and Form factors 

In Bethe Ansatz solvable models, Slavnov's formula 
[TU] is an economical representation of the overlap of an 
eigenstate \{w}) with a general Bethe state \{v}} con- 
structed using the same operators, but for which the set 
of rapidities {v} does not fulfil the Bethe equations. This 
overlap is given as a determinant of an N r by N r matrix, 
which in the problem at hand reads |17j : 

, r 1ir 1V det Nr J({v a },{w b }) U%h,(v b -w a ) 
(W W) = — i=f 7 vrf — ( ' ( 13 ) 

with the matrix elements of J given in Ref . [T7J [H] , 



Jab — 



Vb — Wb 
V a - W b 

-2£ 



E 



l 



(V a — V c )(w b — V c ) 



(14) 



The solution to the inverse problem allows a determi- 
nant representation for the necessary form factors |17j 



(MISSIW) = n 



\T-Q{a)) 



-JL (V a - E a ) Y[(w b - W a ) ]J(v b - Va ) 



(W b — W a 
b>a b<a 



with the matrix elements of T, Q given by 



Tab — 



2 |J(to c - Vb) 

W a — V b 



y 1 -y 1 



Qab(a) = 



(15) 



We explicitly used the fact that both states are solutions 
to the Richardson equations in order to write the matrix 
elements of T in a more compact form than in previous 
publications [T71[T8] . 

The form factors for S~ Sg can be written as sum of 
N r determinants by generalizing the method of Ref. [18] 
for {v} = {w} starting from the double sums in Ref. |17j . 
For a^/3we have 



n 



(M\s-s+\{w}) = 



(W k ~ Wg) 



k^q 



Y[( v a ~ V b ) Y[( w a - Wb) 



b>a 



a>b 



q=l 



(16) 



where, defining A a b = J a bY\(v c — Wb), the matrix ele- 
ments are given by 



<!J ab — A a 



n 



k^b,q 



(w k — W b ) 



9-1 



<J aq 

Jab 



Aaq-1 + 2 



U k -tb+l,q( W k - Wb+l) 

— Aab+1, b 

■OL 

(w q - e fj ){w q - 1 - e a ) 



x Aab+i, b<q-l 

Wb+l — €a 



n 



Y[ (Wk — W q -i) 



Vc - eg 

Wq-l — W q ^ \W C — ep 

(2v a - e a — ep) 



k^q-l 

l/(w a — e a ) 2 , 
A ab , b>q. 



(Va — Ca) 2 (Va — £p) 2 



(17) 



For a = P they are calculated using Hellmann-Feynman 
theorem as explained in the main text. 
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